Algorytm Rungego-Kutty

Algorytm Rungego-Kutty
Metoda Rungego-Kutty

Rodzaj	Metoda numeryczna do iteracyjnego rozwiązywania równań różniczkowych zwyczajnych

Porównanie dokładności metod Rungego-Kutty zastosowanej do rozwiązania równania różniczkowego $y'=\sin(t)^{2}\!\cdot \!y{:}$ czerwony wykres – dokładne rozwiązanie, żółty – metoda RK 4. rzędu, zielony – metoda Heuena, niebieski i różowy – metody Eulera.

Algorytm Rungego-Kutty (metoda Rungego-Kutty) – metoda numeryczna do iteracyjnego rozwiązywania równań różniczkowych zwyczajnych (opisujących procesy fizyczne, chemiczne, biologiczne itp.) Opracowana około 1900 przez niemieckich matematyków: Carla Rungego oraz Martina Kuttę.

Obecnie nazwą metody Rungego-Kutty określa się rodzinę jawnych i niejawnych metod wielokrokowych, jak również pewne ich modyfikacje. Istnieje wiele metod RK, o wielu stopniach, wielu krokach, różnych rzędach, i różniących się między sobą własnościami (jak stabilność, jawność, niejawność, metody osadzone, szybkość działania itp.).

Potocznie metodą Rungego-Kutty określa się najczęściej metodę Rungego-Kutty 4. rzędu ze współczynnikami podanymi poniżej, gdyż jest powszechnie stosowana ze względu na prostotę implementacji, relatywnie proste wzory, dużą szybkość oraz wysoki rząd metody.

Na początku artykułu zostanie omówione podejście do rozwiązywania równań różniczkowych rzędu 1-go, a następnie zostanie pokazane uogólnienie na równania rzędu 2-go i wyższych. Całość dopełnia przykład programu napisany w języku C++, rozwiązujący równanie różniczkowe 2-go rzędu nieliniowe (opis drgań wahadła).

Równanie różniczkowe zwyczajne rzędu 1

Niech dane będzie równanie różniczkowe zwyczajne 1-go rzędu postaci

$y'=f(x,y)$

oraz warunek początkowy

$x_{0},\,y_{0}\equiv y(x_{0}),$

gdzie:

• $x$ – zmienna niezależna,
• $y=y(x)$ – szukana funkcja zmiennej $x,$
• $y'$ – pochodna zmiennej zależnej $y$ względem zmiennej niezależnej $x,$
• $f(x,y)$ – znane wyrażenie algebraiczne, zależne od zmiennych $x$ oraz $y.$

Metoda numeryczna znalezienia rozwiązania polega na obliczeniu dyskretnych wartości $y_{1},y_{2},\dots ,y_{n}$ szukanej funkcji w kolejnych krokach obliczeń. Niech $h$ oznacza przyjęty do obliczeń, niewielki skok zmiennej niezależnej $x,$ tj. zmienna $x$ przyjmuje w kolejnych krokach dyskretne wartości $x_{n},n=0,1,\dots$ zwiększające się o wartość $h.$

Metoda RK 4. rzędu

Iterację zaczynamy zadając $x_{0},\,y_{0}$ oraz funkcję $f(x,y).$ W $n$-tym kroku iteracji, dla $n=0,1,2,\dots ,$ obliczamy następujące wielkości:

(1) $x_{n+1}=x_{n}+h,$

(2) $k_{1}=h\cdot f\left(x_{n},y_{n}\right),$

(3) $k_{2}=h\cdot f{\Big (}x_{n}+{\frac {h}{2}},y_{n}+{\frac {k_{1}}{2}}{\Big )},$

(4) $k_{3}=h\cdot f{\Big (}x_{n}+{\frac {h}{2}},y_{n}+{\frac {k_{2}}{2}}{\Big )},$

(5) $k_{4}=h\cdot f\left(x_{n}+h,y_{n}+k_{3}\right),$

następnie liczymy

(6) $dy_{n}={\frac {1}{6}}(k_{1}+2k_{2}+2k_{3}+k_{4})$

i ostatecznie otrzymamy:

(7) $y_{n+1}=y_{n}+dy_{n}.$

Widać, że $y_{n+1}$ zależy od wartości wcześniej obliczonej $y_{n}$ oraz wielkości kroku $h.$

W ten sposób otrzymuje się (podobnie jak w innych iteracyjnych metodach rozwiązywania równań różniczkowych) kolejne wartości $y_{n},n=1,2,\dots ,$ które przybliżają szukaną funkcję $y.$

Metoda RK 2. rzędu

Istnieje cała klasa metod RK 2. rzędu, z czego z nazwy wyróżnia się przynajmniej dwie metody.

Jawna metoda RK 2, znana jako metoda punktu pośredniego (ang. midpoint method)

$y_{n+1}=y_{n}+h\cdot f\left(x_{n}+{\frac {h}{2}},y_{n}+{\frac {1}{2}}k_{1}\right).$

Oznaczenia takie same, jak wyżej.

Metoda RK 1. rzędu, czyli metoda Eulera

Metoda Eulera jest szczególnym przypadkiem metod Rungego-Kutty (metoda Eulera pojawiła się historycznie najpierw, w związku z tym zachowano tradycyjną nazwę):

$y_{n+1}=y_{n}+k_{1}.$

Oznaczenia takie same jak poprzednio.

Równanie różniczkowe zwyczajne rzędu 2

Jeżeli mamy równanie różniczkowe zwyczajne rzędu 2-go, to zamienia się je na układ dwóch równań rzędu 1-go i stosuje wyżej wymienione metody Rungego-Kutty, iteracyjnie rozwiązując dany układ równań.

Przykład: Równanie drgań wahadła

1. Zamiana równania 2-go rzędu na układ 2 równań 1-go rzędu

Równanie drgań wahadła matematycznego ma postać:

$\theta (t)''+k\cdot \sin \theta (t)=0,$

gdzie:

$\theta (t)$ – kąt odchylenia wahadła od pionu w chwili $t,$

$k=g/\ell$ – stała liczba ($g$ – przyspieszenie ziemskie, $\ell$ – długość nici wahadła).

Podstawiając $\omega =\theta (t)',$ równanie to sprowadza się do układu dwóch równań 1-go rzędu

${\begin{cases}\theta (t)'=\omega ,\\\omega '+k\cdot \sin \theta =0,\end{cases}}$

który można rozwiązać jedną z metod Rungego-Kutty.

2. Obliczenia numeryczne – obliczenie okresu drgań wahadła metodą Eulera

Poniżej podano kod programu w C++ liczącego okres drgań wahadła – zastosowano całkowanie równania ruchu metodą Eulera. Pomimo prostoty metody uzyskuje się bardzo duże dokładności obliczeń, gdy dobierze się odpowiednio mały krok iteracji $dt.$ W programie podanym niżej przejęto $dt=0.000001.$ (Dokładność metody można sprawdzić, porównując wyniki z obliczeniami okresu drgań wahadła, opartymi na całkach eliptycznych pierwszego rodzaju – por. Tabela okresów drgań wahadła)

Program można testować online, korzystając np. z darmowych compilatorów języka C++.

// Metoda Eulera całkowania równania różniczkowego drgań wahadła matematycznego
#include <iostream>
#include <cmath>
using namespace std;

const double g = 9.81, l = 0.25;//9.81;

// Funkcja liczy połowę okresu drgań wahadła matematycznego
double PolOkresDrgan(double theta0, double dt) {
    double theta, theta_2, v, v_2, t, PolOkres;
    // Inicjalizacja zmiennych
    theta = theta0; //Kąt odchylenia od pionu
    v = 0.0;  //prędkość kątowa
    t = 0.0;
    // Obliczania Theta_2, v_2 dla kolejnych kroków czasowych
    while (true) {
        theta_2 = theta + dt * v;
        v_2 = v - dt * (g / l) * sin(theta);
        if(v * v_2 < 0) break; //wahadło osiąga makymalne wychylenie, gdy v zmiania znak - upływa wtedy pół okresu drgań
        // Aktualizacja wartości
        theta = theta_2;
        v = v_2;
        t = t + dt;
    }
    PolOkres = t;
    return PolOkres;
}

int main() {
    double angle, theta0, dt, Okres;
    angle= 179;   //Amplituda drgań w stopniach
    dt = 0.000001; // krok iteracji

    theta0 = angle * M_PI / 180.0;    // Konwersja kąta na radiany
    Okres = 2.0 * PolOkresDrgan(theta0, dt);
    cout << " theta0 = " << angle << " [stopnie]      T = " << Okres <<" [sekundy]"<< endl;

    return 0;
}

Przykład: Obliczenia numeryczne równania rzędu 2-go metodą Rungego-Kutty rzędu 4.

Napisanie programu, rozwiązującego równanie różniczkowe rzędu 2-go metodą RK rzędu 4. wymaga liczenia podwójnej liczby współczynników. Poniżej podano kod programu w języku python; program rozwiązuje równanie różniczkowe wahadła z tłumieniem. Współczynniki metody RK są zdefiniowane w liniach 14 – 30 programu.

Znalezienie rozwiązania analitycznego równania ruchu złożonych układów fizycznych jest w ogólnym przypadku niemożliwe. Jednak metody numeryczne pozwalają efektywnie rozwiązywać te równania ruchu. Np. równanie ruchu wahadła z tłumieniem i z siłą wymuszającą ma postać

${\frac {d^{2}\theta (t)}{dt^{2}}}+2\beta {\frac {d\theta (t)}{dt}}+{\frac {g}{\ell }}\sin \theta (t)=F(t),$

gdzie:

$\beta$ – współczynnik tłumienia,

$F(t)$ – siła wymuszająca, zależna dowolnie od czasu.

Wprowadzając nową zmienną $\omega =\theta (t)'$ powyższe równanie sprowadza się do układu dwóch równań różniczkowych pierwszego rzędu

${\begin{cases}{\frac {d\theta (t)}{dt}}=\omega (t)\\{\frac {d\omega (t)}{dt}}=-2\beta \omega (t)-{\frac {g}{\ell }}\sin \theta (t)+F(t)\end{cases}}$

Następnie układ tych równań rozwiązuje się iteracyjnie, np. metodą Rungego-Kutty, co prowadzi do znalezienia dyskretnych wartości $\theta (t_{n}),\omega (t_{n})$ w zadanym przedziale całkowania równań.

Pokazany tu przykładowy program zakłada zerową siłę wymuszającą, tj. $F(t)=0.$ Uwzględnienie niezerowej siły wymuszającej w kodzie programu jest nieskomplikowane – wystarczy uzupełnić linię 11 kodu, uzupełniając funkcję o dodatkowy parametr – aktualną wartość siły wymuszającej.

Dodatkowo program generuje cztery wykresy zależności $T(\theta _{0},c),$ gdzie $c$ – współczynnik tłumienia, korzystając online z biblioteki mathplotlib.pyplot.

Program można testować, korzystając np. z darmowego notatnika colab google online.

import math
import matplotlib.pyplot as plt

g = 9.81  # Przyspieszenie ziemskie
L = 9.81   # Długość wahadła

T_0 = 2 * math.pi * math.sqrt(L / g)

def f(theta, omega, c):
    # Równania ruchu dla wahadła z tłumieniem
    omega_dot = (-c * omega - (g / L) * math.sin(theta))
    return omega_dot

def runge_kutta(theta, omega, h, c):
    k1_theta = h * omega
    k1_omega = h * f(theta, omega, c)

    k2_theta = h * (omega + k1_omega / 2)
    k2_omega = h * f(theta + k1_theta / 2, omega + k1_omega / 2, c)

    k3_theta = h * (omega + k2_omega / 2)
    k3_omega = h * f(theta + k2_theta / 2, omega + k2_omega / 2, c)

    k4_theta = h * (omega + k3_omega)
    k4_omega = h * f(theta + k3_theta, omega + k3_omega, c)

    theta += (k1_theta + 2 * k2_theta + 2 * k3_theta + k4_theta) / 6
    omega += (k1_omega + 2 * k2_omega + 2 * k3_omega + k4_omega) / 6

    return theta, omega

def main():
    # DANE:
    theta_0 = 170  # Początkowy kąt wychylenia [w stopniach]

    h = 0.00001          # Krok czasu
    iterations = 4000000  # Liczba iteracji

    c_values = [0, 0.2, 0.5, 1.5]  # Wartości współczynnika tłumienia c
    print(c_values)

    for c in c_values:
        omega = 0
        theta = theta_0 / 180 * math.pi
        time_values = []
        theta_values = []

        print("Tłumienie c =", c)
        time_2 = 0
        for i in range(iterations):
            time = i * h
            time_values.append(time / T_0)
            theta_values.append(theta)

            omega_2 = omega
            theta, omega = runge_kutta(theta, omega, h, c)
            if (omega * omega_2 < 0) or (omega_2 == 0):
                print("T/2/T_0=", (time - time_2) / T_0)
                time_2 = time

        plt.plot(time_values, theta_values)

    # WYKRESY - RYSOWANIE:
    plt.xlabel(r'Czas unormowany $t/T_0$')
    plt.ylabel(r'$\theta(t)$ [stopnie]')
    plt.title('Wykres kąta wychylenia wahadła z tłumieniem od czasu')
    plt.grid(True)
    plt.legend(['c = 0', 'c = 0.2', 'c = 0.5', 'c = 1.5'])

    plt.show()

if __name__ == "__main__":
    main()

Równanie różniczkowe zwyczajne rzędu n

Jeżeli mamy równanie różniczkowe zwyczajne rzędu $n,$ to w szczególnych wypadkach można obniżyć jego rząd. Zagadnienie obniżania stopnia równań różniczkowych omawia np. Np. gdy dane jest równanie różniczkowe w postaci jawnej

$F\left(x,y,y',y'',\ \dots ,\ y^{(n-1)}\right)=y^{(n)},$

to zamienia się je na układ $n$ równań rzędu 1 poprzez wprowadzenie dodatkowych $n-1$ zmiennych $u_{i}=y^{(i-1)},\quad i=2,3,\dots ,n,$

tj. pochodne wyższych rzędów niż 1 traktuje się jako nowe zmienne; przyjmując dodatkowo oznaczenie $u_{1}\equiv y$ dla jednolitości zapisu otrzymamy układ n równań różniczkowych pierwszego rzędu w postaci:

${\begin{array}{rcl}u_{1}'&=&u_{2}\\u_{2}'&=&u_{3}\\&\vdots &\\u_{n-1}'&=&u_{n}\\u_{n}'&=&F(x,u_{1},\dots ,u_{n}).\end{array}}$

lub bardziej zwięźle w notacji wektorowej:

$\mathbf {u} '=\mathbf {F} (x,\mathbf {u} )$

gdzie:

$\mathbf {u} =(u_{1},\dots ,u_{n}),\quad \mathbf {F} (x,u_{1},\dots ,u_{n})=[u_{2},u_{3},\dots ,u_{n},F(x,u_{1},\dots ,u_{n})]$

Otrzymany układ równań rozwiązuje się stosując jedną z wyżej wymienionych metod Rungego-Kutty.

Zobacz też

Całkowanie równań różniczkowych:

• metoda Eulera
• metoda Eulera-Maruyamy
• równanie sztywne - problem stabilności rozwiązań numerycznych

Inne:

• całkowanie numeryczne
• metody numeryczne
• metoda Monte Carlo
• metoda Monte-Carlo Metropolisa-Hastingsa

Bibliografia

• Z. Fortuna, B. Macukow, J. Wąsowski, Metody numeryczne Podręczniki akademickie Elektronika, informatyka, telekomunikacja, Warszawa: Wydawnictwa Naukowo-Techniczne, 1982, s. 285–312.
• R.S. Guter, A.R. Janpolski: Równania różniczkowe. Warszawa: Wydawnictwo Naukowe PWN, 1980.
• D. Potter, Metody obliczeniowe fizyki, fizyka komputerowa, Warszawa: PWN, 1982, s. 19–43.
• J. Szmelter, Metody komputerowe w mechanice, Warszawa: PWN, 1980, s. 150–157.

Linki zewnętrzne

• MyPhysicsLab – Runge-Kutta Algorithm – omówienie (w języku angielskim)
• Numerical Recipes in C – Runge-Kutta Method – algorytm w języku C